Generalized Spatial Regression with Partial Differential Equation Regularization

Simone Panzeri @ fdaPDE Team, MOX, Department of Mathematics, Politecnico di Milano, Italy

The fdaPDE C++/R library (Palummo et al., 2025) provides advanced statistical methods for physics-informed spatial and functional data analysis, featuring regularization terms based on partial differential equations; for a review of this class of methods, refer to Sangalli (2021). The library can handle different types of data that may evolve over time on non-standard spatial domains, including linear networks, irregular planar regions, curved surfaces, and volumes. Among the broad range of modeling capabilities offered by fdaPDE, we focus here on the generalized spatial regression method.

Generalized spatial regression method

Let \(\mathcal{D} \subset \mathbb{R}^d\), with \(d \geq 1\), be a bounded spatial domain of interest in which \(n\) fixed measurement stations are located. Here, we consider \(d = 2\), but the proposed method can handle multidimensional spatial domains with complex shapes, including two-dimensional curved surfaces (Ettinger et al., 2016) and non-convex three-dimensional volumes (Arnone et al., 2023). At the location \(\mathbf{p}_i = (p_{1i}, p_{2i})^\top \in \mathcal{D}\) of each station, we observe a realization \(y_i \in \mathbb{R}\) of the response variable \(Y_i\) under study, along with a set of covariates \(\mathbf{x}_i = \left( x_{i,1},\ldots,x_{i,q} \right)^\top \in \mathbb{R}^q\), if available. Specifically, we assume that \(Y_i\) follows a distribution from the exponential family, with mean \(\mu_i\) and a common scale parameter \(\phi\), for \(i = 1, \ldots, n\). In the context of spatial regression, we present a generalized additive model with partial differential equation regularization, which we formulate according to Wilhelm & Sangalli (2016), as follows: \[g(\mu_i) = \theta_i = \mathbf{x}_i^\top \boldsymbol{\beta} + f(\mathbf{p}_i) \,, \qquad i = 1, \ldots, n,\] where:

To estimate the unknown parameters, we refer to the basic form of the regularized methods, as reviewed in Sangalli (2021) and originally introduced in Sangalli et al. (2013). These works propose to maximize the following penalized functional based on the log-likelihood \(\ell(\cdot ; \boldsymbol{\theta})\) for the considered distribution: \[\sum_{i = 1}^n \ell\left( y_i; \theta_i \right) - R(f; \lambda)\,,\] where \(R(f; \lambda)\) is the regularization term that depends on the smoothing parameter \(\lambda \in \mathbb{R}^+\). The regularization term is defined using Partial Differential Equations (PDEs), with differential operators encoding the smoothness of \(f\) over \(\mathcal{D}\). This term may also incorporate problem-specific information into the modeling framework, whenever available. The functional above embodies the trade-off between data fidelity and model fidelity through the least-squares term and the misfit with respect to the PDE, respectively. The balance between these two contributions is governed by the smoothing parameter. When the considered distribution is Gaussian, the optimization of the functional above coincides with the one that follows from the basic form of spatial regression with PDE regularization.

In its basic form, the regularization term is given by: \[R(f; \lambda) = \lambda \int_\mathcal{D} \left( \Delta f (\mathbf{p}) \right)^2 \, \text{d}\mathbf{p}\,,\] where the Laplacian operator is defined as \[\Delta f = \frac{\partial^2 f}{\partial p_1^2} + \frac{\partial^2 f}{\partial p_2^2}\] and is a measure of the local curvature of the field \(f\). More complex roughness penalties may be considered in the regularization term, with the PDE modeling the spatial variation of the unknown spatial field \(f\), thereby accounting for spatial non-stationarities and anisotropies, as discussed in Wilhelm & Sangalli (2016).

The method for Generalized Spatial Regression with Partial Differential Equation regularization presented in this vignette, hereafter referred to as GSR-PDE, is implemented within the fdaPDE C++/R library (Palummo et al., 2025). We load the library in the working environment as follows:

# Import the fdaPDE library
library(fdaPDE2)           # v. 2.0 (2025)
rm(list = ls())

# Load additional libraries and helper functions for plotting
source("../utils/graphics.R")

Application to fire count data in Southern Italy

As a benchmark application in environmental sciences, we consider fire count data recorded in Southern Italy, during 2021. In addition to fire counts per administrative division (province), the dataset includes several areal covariates, such as area, elevation, average temperature, solar radiation, total precipitation, population, land use, and land cover, as of 2021. The dataset is freely and publicly available by different data sources, including the Fire Information for Resource Management System (FIRMS) of the National Aeronautics and Space Administration (NASA), the CORINE Land Cover (CLC) from the European Environment Agency (EEA), the Digital Elevation Model (DEM) from the National Institute of Geophysics and Volcanology (INGV), and the Copernicus ERA5 from the European Centre for Medium-Range Weather Forecasts (ECMWF) as part of the Copernicus Climate Change Service (C3S). Due to ongoing climate change and global warming, large-scale wildfires are occurring with increasing frequency worldwide, especially in southern Europe. This results in severe repercussions and damage to ecosystems, infrastructures, material assets, and socio-economic activities. Therefore, an in-depth statistical analysis of this phenomenon is becoming crucial for driving the implementation of new policies aimed at preserving both artificial and natural environments.

Spatial domain

We import the preprocessed data observed in Southern Italy during 2021 from the shapefile ../data/GSRPDE_2D/GSRPDE_2D_data.shp as a sf (Pebesma, 2025) object.

## [DATA]
# Load the data
data_sf <- st_read(dsn = "../data/GSRPDE_2D/GSRPDE_2D_data.shp")
#> Reading layer `GSRPDE_2D_data' from data source 
#>   `/home/user/IWSM25_short_course/data/GSRPDE_2D/GSRPDE_2D_data.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 33 features and 17 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 12.44208 ymin: 36.64986 xmax: 18.52069 ymax: 42.89375
#> Geodetic CRS:  WGS 84
head(data_sf)
#> Simple feature collection with 6 features and 17 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 12.89095 ymin: 37.09486 xmax: 18.09681 ymax: 41.48737
#> Geodetic CRS:  WGS 84
#>   COUNTRY   REGION              PROVINCE TERRITORY     AREA MIN_ELEV AVG_ELEV
#> 1   Italy  Sicilia             Agrigento    Island 3025.283   0.5303 310.0765
#> 2   Italy Campania              Avellino     South 2803.711  92.7902 610.1539
#> 3   Italy   Puglia                  Bari     South 3797.713   1.0050 302.0524
#> 4   Italy   Puglia Barletta-Andria-Trani     South 1524.347   0.4824 217.2597
#> 5   Italy Campania             Benevento     South 2045.429  33.2261 479.9234
#> 6   Italy   Puglia              Brindisi     South 1822.502   0.4948 109.3009
#>    MAX_ELEV ARTIFIC AGRICUL FOREST WETLANDS  WATER FIRE_COUNT     POP REF_AREA
#> 1 1414.2872  0.0360  0.7901 0.1641   0.0002 0.0096        306  413177  3052.82
#> 2 1695.8716  0.0445  0.6325 0.3213   0.0002 0.0015         16  398932  2805.96
#> 3  659.2282  0.0633  0.8279 0.1074   0.0000 0.0014         75 1225048  3862.66
#> 4  657.4872  0.0448  0.8197 0.1005   0.0303 0.0047         26  379509  1542.99
#> 5 1719.1609  0.0371  0.7076 0.2527   0.0000 0.0026         17  263125  2080.37
#> 6  394.1095  0.0658  0.9121 0.0153   0.0018 0.0050         24  379522  1861.33
#>   POP_DENS                       geometry
#> 1 135.3427 MULTIPOLYGON (((12.89625 37...
#> 2 142.1731 MULTIPOLYGON (((14.60641 40...
#> 3 317.1514 MULTIPOLYGON (((17.35441 40...
#> 4 245.9569 MULTIPOLYGON (((16.20329 40...
#> 5 126.4799 MULTIPOLYGON (((14.5728 41....
#> 6 203.8983 MULTIPOLYGON (((18.09681 40...

# Number of data (i.e., of provinces)
n <- nrow(data_sf)

Before entering the details about the data, we visualize the boundaries of each province within the spatial domain of interest in Southern Italy. These boundaries are already embedded in the sf object loaded above. To enable interactive visualization, we use the mapview (Appelhans et al., 2023) R package.

# Provinces
provinces_sfc <- st_geometry(obj = data_sf)
provinces_sf <- st_as_sf(x = provinces_sfc)
provinces_sf
#> Simple feature collection with 33 features and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 12.44208 ymin: 36.64986 xmax: 18.52069 ymax: 42.89375
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                                 x
#> 1  MULTIPOLYGON (((12.89625 37...
#> 2  MULTIPOLYGON (((14.60641 40...
#> 3  MULTIPOLYGON (((17.35441 40...
#> 4  MULTIPOLYGON (((16.20329 40...
#> 5  MULTIPOLYGON (((14.5728 41....
#> 6  MULTIPOLYGON (((18.09681 40...
#> 7  MULTIPOLYGON (((14.44479 37...
#> 8  MULTIPOLYGON (((14.50486 41...
#> 9  MULTIPOLYGON (((14.03236 40...
#> 10 MULTIPOLYGON (((14.79454 37...

# Boundaries of each province
boundary_sf <- st_boundary(x = provinces_sf)
boundary_sf
#> Simple feature collection with 33 features and 0 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 12.44208 ymin: 36.64986 xmax: 18.52069 ymax: 42.89375
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                                 x
#> 1  MULTILINESTRING ((12.89625 ...
#> 2  MULTILINESTRING ((14.60641 ...
#> 3  MULTILINESTRING ((17.35441 ...
#> 4  MULTILINESTRING ((16.20329 ...
#> 5  MULTILINESTRING ((14.5728 4...
#> 6  MULTILINESTRING ((18.09681 ...
#> 7  MULTILINESTRING ((14.44479 ...
#> 8  MULTILINESTRING ((14.50486 ...
#> 9  MULTILINESTRING ((14.03236 ...
#> 10 MULTILINESTRING ((14.79454 ...

# Interactive plot
mapview(provinces_sf, col.regions = "gray75", lwd = 1.5,
        legend = FALSE, layer.name = "provinces")

Figure 1: Boundaries of each of the 33 provinces in Southern Italy.

Now we build a regular mesh of the spatial domain above. The mesh provides a discretization of the spatial domain, which is required to solve the estimation problem using the finite element method. To do so, we directly load mesh nodes and triangles from ../data/GSRPDE_2D/GSRPDE_2D_mesh_nodes.txt and ../data/GSRPDE_2D/GSRPDE_2D_mesh_triangles.txt files. These files store a matrix of size #nodes-by-2 containing the \(x\) and \(y\) coordinates of the mesh nodes, and a matrix of size #triangles-by-3, where the \(i\)th row contains the indices of the mesh nodes forming the \(i\)th triangle. These matrices are obtained using the femR (Clemente et al., 2025) and RTriangle (Shewchuk & Sterratt, 2025) R packages; further details are omitted here. Other software can be used to get the triangulation; for example, it can be constructed through the fmesher (Lindgren, 2025) R package.

## [MESH]
# Load the mesh nodes
mesh_nodes <- read.table(file = "../data/GSRPDE_2D/GSRPDE_2D_mesh_nodes.txt")
head(mesh_nodes)
#>         V1       V2
#> 1 17.02403 39.48284
#> 2 16.97541 39.49125
#> 3 16.86791 39.53792
#> 4 16.83847 39.55597
#> 5 16.81708 39.58875
#> 6 16.78153 39.61458

# Load the nodes markers
mesh_nodesmarkers <- read.table(file = "../data/GSRPDE_2D/GSRPDE_2D_mesh_nodesmarkers.txt")
head(mesh_nodesmarkers)
#>   V1
#> 1  1
#> 2  1
#> 3  1
#> 4  1
#> 5  1
#> 6  1

# Load the mesh triangles
mesh_triangles <- read.table(file = "../data/GSRPDE_2D/GSRPDE_2D_mesh_triangles.txt")
head(mesh_triangles)
#>    V1   V2   V3
#> 1 934  283  949
#> 2 424 1068 1582
#> 3 283  284  949
#> 4 539 1296 1095
#> 5 273  271  272
#> 6  79   80 1082

# Set up the triangulation for fdaPDE
mesh <- triangulation(cells = mesh_triangles, nodes = mesh_nodes,
                      boundary = mesh_nodesmarkers)

# Nodes coordinates
head(mesh$nodes)
#>          [,1]     [,2]
#> [1,] 17.02403 39.48284
#> [2,] 16.97541 39.49125
#> [3,] 16.86791 39.53792
#> [4,] 16.83847 39.55597
#> [5,] 16.81708 39.58875
#> [6,] 16.78153 39.61458

# Number of nodes
mesh$n_nodes
#> [1] 2563

# Edges
head(mesh$nodes)
#>          [,1]     [,2]
#> [1,] 17.02403 39.48284
#> [2,] 16.97541 39.49125
#> [3,] 16.86791 39.53792
#> [4,] 16.83847 39.55597
#> [5,] 16.81708 39.58875
#> [6,] 16.78153 39.61458

# Number of edges
mesh$n_edges
#> [1] 7075

# Triangles
head(mesh$cells)
#>      [,1] [,2] [,3]
#> [1,]  934  283  949
#> [2,]  424 1068 1582
#> [3,]  283  284  949
#> [4,]  539 1296 1095
#> [5,]  273  271  272
#> [6,]   79   80 1082

# Number of triangles
mesh$n_cells
#> [1] 4514

# Bounding Box
mesh$bbox
#>          [,1]     [,2]
#> [1,] 12.44208 36.64986
#> [2,] 18.52069 42.89375

We visualize the resulting mesh of the spatial domain of interest.

# Interactive plot
mapview(boundary_sf,
        col.regions = "gray25", alpha.regions = 0.25, col = "black", lwd = 1.5,
        legend = FALSE, layer.name = "domain") +
  mapview(mesh, crs = 4326, col.regions = "transparent", lwd = 1.25,
          legend = FALSE, layer.name = "mesh")

Figure 2: Regular mesh of the spatial domain of interest of Southern Italy with \(2\,568\) nodes and \(4\,523\) triangles. Note that the mesh conforms to the boundaries of the provinces.

We use the triangulation just created to define the spatial support of a geoframe object. This object will host layers containing data observed over the spatial support.

# Create the geoframe
italy <- geoframe(domain = mesh)
italy
#> Geoframe with 0 layers
#> Bounding box:    xmin: 12.44208 ymin: 36.64986 xmax: 18.52069 ymax: 42.89375
#> Number of nodes: 2563
#> Number of cells: 4514

Data

We briefly introduce the data under study. We add a data layer to the geoframe object defined above, thus obtaining a format compatible with the implementation of the proposed regression method.

# Add layer with data to the geoframe object (directly from the shapefile)
italy$load_shp(layer = "fires", filename = "../data/GSRPDE_2D/GSRPDE_2D_data.shp")
italy
#> Geoframe with 1 layers
#> Bounding box:    xmin: 12.44208 ymin: 36.64986 xmax: 18.52069 ymax: 42.89375
#> Number of nodes: 2563
#> Number of cells: 4514
#> 
#> Layer: fires
#> Type:  AREAL
#> Dims: 33, 17
#> First 6 data rows:
#> COUNTRY REGION   PROVINCE              TERRITORY AREA    MIN_ELEV AVG_ELEV MAX_ELEV ARTIFIC AGRICUL FOREST  WETLANDS WATER   FIRE_COUNT POP     REF_AREA POP_DENS 
#> <chr>   <chr>    <chr>                 <chr>     <flt64> <flt64>  <flt64>  <flt64>  <flt64> <flt64> <flt64> <flt64>  <flt64> <flt64>    <flt64> <flt64>  <flt64>  
#> Italy   Sicilia  Agrigento             Island    3025.28 0.5303   310.076  1414.29  0.036   0.7901  0.1641  2e-04    0.0096  306        413177  3052.82  135.343  
#> Italy   Campania Avellino              South     2803.71 92.7902  610.154  1695.87  0.0445  0.6325  0.3213  2e-04    0.0015  16         398932  2805.96  142.173  
#> Italy   Puglia   Bari                  South     3797.71 1.005    302.052  659.228  0.0633  0.8279  0.1074  0        0.0014  75         1225048 3862.66  317.151  
#> Italy   Puglia   Barletta-Andria-Trani South     1524.35 0.4824   217.26   657.487  0.0448  0.8197  0.1005  0.0303   0.0047  26         379509  1542.99  245.957  
#> Italy   Campania Benevento             South     2045.43 33.2261  479.923  1719.16  0.0371  0.7076  0.2527  0        0.0026  17         263125  2080.37  126.48   
#> Italy   Puglia   Brindisi              South     1822.5  0.4948   109.301  394.11   0.0658  0.9121  0.0153  0.0018   0.005   24         379522  1861.33  203.898

The variables observed at the province level are given by:

In addition, “geometry” contains the multipolygons defining the provinces under study.

The values of most of these variables, recorded by province in 2021, are shown in the interactive plot below, created using mapview.

# Color palettes
color_palette_fire_counts <- c("#FFFFB2", "#FECC5C", "#FD8D3C", "#F03B20", "#BD0026")
color_palette_population <- c("#DEEBF7", "#9ECAE1", "#6BAED6", "#3182BD", "#08519C")
color_palette_forest <- c("#E5F5E0", "#A1D99B", "#74C476", "#31A354", "#006D2C")
color_palette_elevation <- c("#006400", "#A2CD5A", "#F5DEB3", "#8B864E", "#8B2323")

# Interactive plot
mapview(italy[["fires"]], crs = 4326, varnames = c("FIRE_COUNT", "POP", "FOREST", "AVG_ELEV"),
        color_palettes = list(color_palette_fire_counts,
                              color_palette_population,
                              color_palette_forest,
                              color_palette_elevation))

Figure 3: Fire counts observed at 33 provinces in Southern Italy during 2021 (top left); population measured at the province level as of 2021 (top right); forest and seminatural areas (in \(\text{km}^2\)) measured at the province level as of 2021 (bottom left); average elevation (in \(\text{m}\)) measured at the province level as of 2021 (bottom right).

First, we focus on the FIRE_COUNT variable. The non-standard shape of the spatial domain strongly influences the quantity under study. For example, fire counts recorded in the two provinces facing the Strait of Messina, namely the provinces of Messina in Sicily and Reggio-Calabria in Calabria, are not highly correlated. This spatial variation suggests that fire counts depends on factors beyond purely environmental ones (e.g., temperature, elevation, land cover), which typically vary smoothly across space. Specifically, differences in fire counts may be partly due to the varying number of arson incidents reported in 2021 in the two provinces, highlighting the role of unobserved cultural and social phenomena, which may not be fully captured by methods that rely solely on Euclidean distances. Differently, GSR-PDE naturally accounts for the geometry of the domain, which may feature disjoint regions with irregular boundaries and concavities, as in the case under consideration.

GSR-PDE model fitting for Poisson response (without covariates)

Now we are ready to perform the GSR-PDE spatial regression method using the gsr function. This function offers several options for solving the regression problem, including the possibility to specify covariates, PDE parameters, and boundary conditions, whenever available. The function supports several distributions within the exponential family for generalized linear models and can handle areal data, as in the case under study. In addition, the function implements different algorithms for degrees of freedom computation and criteria for optimal smoothing parameter selection. For further information, see the documentation in the Help by typing ?gsr.

Poisson regression and intensity estimation

We now show the pointwise spatial prediction. To evaluate the regression model fit over a fine grid and enable interactive visualization of the estimate, we internally create a new raster object. We then compute the fitted values over the grid using the $eval method. We plot the resulting estimate with mapview: this is just one possibility; various other plotting options are available depending on the specific purpose.

The regression model fit over the fine grid is displayed below (left). In addition, we load and show the point pattern associated with the locations of fires occurred in 2021 (right).

# Load locations of fires
locations <- st_read(dsn = "../data/GSRPDE_2D/GSRPDE_2D_locations.shx")
#> Reading layer `GSRPDE_2D_locations' from data source 
#>   `/home/user/IWSM25_short_course/data/GSRPDE_2D/GSRPDE_2D_locations.shx' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 4510 features and 5 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 12.5088 ymin: 36.7219 xmax: 18.4792 ymax: 42.8807
#> Geodetic CRS:  WGS 84

# Interactive plot
map1 <- mapview(f_grid, crs = 4326, col.regions = color_palette_fire_counts,
                alpha.regions = 0.75, na.color = "transparent",
                layer.name = "intensity") +
  mapview(boundary_sf, color = "gray25", lwd = 1.5,
          legend = FALSE, layer.name = "provinces")

map2 <- mapview(locations, legend = FALSE, col.regions = "black",
                alpha.regions = 0.75, stroke = FALSE, cex = 2,
                layer.name = "locations") +
  mapview(boundary_sf, color = "gray25", lwd = 1.5,
          legend = FALSE, layer.name = "provinces")

sync(map1, map2)

Figure 6: Estimated spatial field of the fire counts provided by the isotropic GSR-PDE (without covariates) with \(\lambda\) selected via GCV minimization using grid search (left); locations of fires occurred in Southern Italy in 2021 (right).

In particular, by focusing on the spatial point pattern rather than aggregated counts, the Poisson regression and the inhomogeneous Poisson point process share equivalent modeling structures, enabling a direct comparison of their pointwise intensity estimates. For density and intensity estimation with PDE regularization, see Ferraccioli et al. (2021) and Begu et al. (2024).

GSR-PDE model fitting for Poisson response (future developments)

The phenomenon under study is influenced by several factors, including the average elevation, the population, the type of territory, as well as the amount of artificial areas, agricultural areas, forest and seminatural areas, wetlands and water bodies, all measured at the province level as of 2021. The proposed GSR-PDE method can accommodate these covariates into the modeling framework. We consider the three covariates shown above – population, average elevation (in \(\text{m}\)), and forest and seminatural areas (in \(\text{km}^2\)) – by province, as of 2021. These covariates could be included into the modeling framework of GSR-PDE to further explain the fire count data in Southern Italy.

Moreover, as a future development, this model fitting could leverage problem-specific information derived from the physics of underlying factors influencing the phenomenon considered here. In particular, the regularization term could incorporate a partial differential equation accounting for the presence of wind across the provinces.

References

Appelhans, T., Detsch, F., Reudenbach, C., & Woellauer, S. (2023). Mapview: Interactive viewing of spatial data in R. url: https://CRAN.R-project.org/package=mapview.
Arnone, E., Negri, L., Panzica, F., & Sangalli, L. M. (2023). Analyzing data in complicated 3D domains: Smoothing, semiparametric regression, and functional principal component analysis. Biometrics, 79(4), 3510–3521. https://doi.org/10.1111/biom.13845
Begu, B., Panzeri, S., Arnone, E., Carey, M., & Sangalli, L. M. (2024). A nonparametric penalized likelihood approach to density estimation of space-time point patterns. Spatial Statistics. https://doi.org/10.1016/j.spasta.2024.100824
Clemente, A., Palummo, A., Arnone, E., Formaggia, L., & Sangalli, L. M. (2025). femR: Bridging physics and statistics in R. url: https://github.com/fdaPDE/femR.
Ettinger, B., Perotto, S., & Sangalli, L. M. (2016). Spatial regression models over two-dimensional manifolds. Biometrika, 103(1), 71–88.
Ferraccioli, F., Arnone, E., Finos, L., Ramsay, J. O., & Sangalli, L. M. (2021). Nonparametric density estimation over complicated domains. Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(2), 346–368. https://doi.org/10.1111/rssb.12415
Lindgren, F. (2025). Fmesher: Triangle meshes and related geometry tools. url: https://CRAN.R-project.org/package=fmesher.
Palummo, A., Clemente, A., Sangalli, E., Arnone, Formaggia, L., & Maria, L. (2025). fdaPDE: Physics-informed statistical learning. url: https://CRAN.R-project.org/package=fdaPDE.
Pebesma, E. (2025). Sf: Simple Features for R. url: https://CRAN.R-project.org/package=sf.
Sangalli, L. M. (2021). Spatial regression with partial differential equation regularisation. International Statistical Review, 89(3), 505–531. https://doi.org/10.1111/insr.12444
Sangalli, L. M., Ramsay, J. O., & Ramsay, T. O. (2013). Spatial spline regression models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 75(4), 681–703. https://doi.org/10.1111/rssb.12009
Shewchuk, J. R., & Sterratt, D. C. (2025). RTriangle - a 2D quality mesh generator and Delaunay triangulator. url: https://CRAN.R-project.org/package=RTriangle.
Wilhelm, M., & Sangalli, L. M. (2016). Generalized spatial regression with differential regularization. Journal of Statistical Computation and Simulation, 86(13), 2497–2518. https://doi.org/10.1080/00949655.2016.1182532